Linear-scaling quantum Monte Carlo technique with non-orthogonal localized orbitals

We have reformulated the quantum Monte Carlo (QMC) technique so that a large part of the calculation scales linearly with the number of atoms. The reformulation is related to a recent alternative proposal for achieving linearscaling QMC, based on maximally localized Wannier orbitals (MLWO), but has the advantage of greater simplicity. The technique we propose draws on methods recently developed for linear-scaling density functional theory. We report tests of the new technique on the insulator MgO, and show that its linear-scaling performance is somewhat better than that achieved by the MLWO approach. Implications for the application of QMC to large complex systems are pointed out. The quantum Monte Carlo (QMC) technique [1] is becoming ever more important in the study of condensed matter, with recent applications including the reconstruction of semiconductor surfaces [2], the energetics of point defects in insulators [3], optical excitations in nanostructures [4], and the energetics of organic molecules [5]. Although its demands on computer power are much greater than those of widely used techniques such as density functional theory (DFT), its accuracy is also much greater for most systems. With QMC now being applied to large complex systems containing hundreds of atoms, a major issue is the scaling of the required computer effort with system size. In other electronic-structure techniques, including DFT, the locality of quantum coherence [6] suggests that it should generally be possible to achieve linear-scaling,or O(N) operation, in which the computer effort is proportional to the number of atoms N . Very recently, a procedure has been suggested [7] for achieving at least partial linear scaling for QMC, based on the idea of ‘maximally localized Wannier functions’ [8]. The purpose of this letter is to propose and test a simpler alternative method, which appears to have important advantages. 0953-8984/04/250305+07$30.00 © 2004 IOP Publishing Ltd Printed in the UK L305 L306 Letter to the Editor The O(N) techniques that have been developed for tight binding (TB) [9],DFT [10, 11] and Hartree–Fock [12] calculations all depend ultimately on the fact that the density matrix ρ(r, r′) associated with the single-electron orbitals decays to zero as |r − r′| → ∞, and the manner of this decay has been extensively studied ([13] and references therein). Briefly, the decay is algebraic for metals and exponential for insulators, with the decay rate increasing with band gap, so that there is more to be gained from O(N) techniques for wide-gap insulators. Equivalently, the extended orbitals used in most conventional techniques can be linearly combined to form localized orbitals, which again decay exponentially in insulators. (For a review of early localized-orbital methods in quantum chemistry, see [14].) Orthogonal Wannier functions are one form of localized orbitals, but it has long been recognized that stronger localization can be achieved by going to non-orthogonal orbitals, as is done in many existing O(N) TB, DFT and quantum-chemistry techniques. In QMC, the trial many-body wavefunction T(r1, . . . , rN ) consists of a Slater determinant D of single-electron orbitals ψn(ri) multiplied by a parameterized Jastrow correlation factor J (r1, . . . , rN ). (In the pseudopotential-based QMC of interest here, the ψn(ri) are commonly taken from a plane-wave pseudopotential DFT calculation.) In variational Monte Carlo (VMC), J is ‘optimized’ by varying its parameters so as to reduce the variance of the ‘local energy’ −1 T (Ĥ T), where Ĥ is the many-electron Hamiltonian. Since VMC by itself is not usually accurate enough, the optimized T produced by VMC is used in diffusion Monte Carlo (DMC), which achieves the exact ground state within the fixed nodal structure imposed by the Slater determinant D. In conventional DMC, a large fraction of the computer time goes into evaluating the single-electron orbitals for all the electron positions ri in each of the replicas (QMC ‘walkers’). For this part of the calculations, the number of computer operations required to perform each QMC step scales at least as N2 (M orbitals ψn(ri) for N positions ri , with M proportional to N), and the scaling deteriorates to N3 if (as is often done) a plane-wave basis set is used to represent the ψn(ri). The memory requirement scales as N2, and this also places important limits on the size of system that can be treated. As pointed out by Williamson et al [7], the scaling for the evaluation of the ψn(ri) can be reduced from N3 or N2 to N if the Slater determinant D is re-expressed in terms of localized orbitals, which themselves are represented in terms of localized basis functions. The key point here is that a determinant is changed by at most a constant overall factor if arbitrary linear combinations of its rows or columns are made. More precisely, if we construct orbitals φm(r) as linear combinations